
/***************************************************************************
                                                                           *
Copyright 2013 CertiVox UK Ltd.                                           *
                                                                           *
This file is part of CertiVox MIRACL Crypto SDK.                           *
                                                                           *
The CertiVox MIRACL Crypto SDK provides developers with an                 *
extensive and efficient set of cryptographic functions.                    *
For further information about its features and functionalities please      *
refer to http://www.certivox.com                                           *
                                                                           *
* The CertiVox MIRACL Crypto SDK is free software: you can                 *
  redistribute it and/or modify it under the terms of the                  *
  GNU Affero General Public License as published by the                    *
  Free Software Foundation, either version 3 of the License,               *
  or (at your option) any later version.                                   *
                                                                           *
* The CertiVox MIRACL Crypto SDK is distributed in the hope                *
  that it will be useful, but WITHOUT ANY WARRANTY; without even the       *
  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
  See the GNU Affero General Public License for more details.              *
                                                                           *
* You should have received a copy of the GNU Affero General Public         *
  License along with CertiVox MIRACL Crypto SDK.                           *
  If not, see <http://www.gnu.org/licenses/>.                              *
                                                                           *
You can be released from the requirements of the license by purchasing     *
a commercial license. Buying such a license is mandatory as soon as you    *
develop commercial activities involving the CertiVox MIRACL Crypto SDK     *
without disclosing the source code of your own applications, or shipping   *
the CertiVox MIRACL Crypto SDK with a closed source product.               *
                                                                           *
***************************************************************************/
/*
 *   MIRACL Montgomery's method for modular arithmetic without division.
 *   mrmonty.c
 *
 *   Programs to implement Montgomery's method
 *   See "Modular Multiplication Without Trial Division", Math. Comp.
 *   Vol 44, Number 170, April 1985, Pages 519-521
 *   NOTE - there is an important correction to this paper mentioned as a
 *   footnote in  "Speeding the Pollard and Elliptic Curve Methods",
 *   Math. Comput., Vol. 48, January 1987, 243-264
 *
 *   The advantage of this approach is that no division required in order
 *   to compute a modular reduction - useful if division is slow
 *   e.g. on a SPARC processor, or a DSP.
 *
 *   The disadvantage is that numbers must first be converted to an internal
 *   "n-residue" form.
 *
 */

#include <stdlib.h>

#include "miracl.h"

#ifdef MR_FP
#include <math.h>
#endif

#ifdef MR_WIN64
#include <intrin.h>
#endif

#ifdef MR_COUNT_OPS
extern int fpc, fpa;
#endif

#ifdef MR_CELL
extern void mod256(_MIPD_ big, big);
#endif

void kill_monty(_MIPDO_)
{
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif
    zero(mr_mip->modulus);
#ifdef MR_KCM
    zero(mr_mip->big_ndash);
#endif
}

mr_small prepare_monty(_MIPD_ big n)
{ /* prepare Montgomery modulus */
#ifdef MR_KCM
    int nl;
#endif
#ifdef MR_PENTIUM
    mr_small ndash;
    mr_small base;
    mr_small magic   = 13835058055282163712.0;
    int      control = 0x1FFF;
#endif
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif
    if (mr_mip->ERNUM) return (mr_small)0;
    /* Is it set-up already? */
    if (size(mr_mip->modulus) != 0)
        if (mr_compare(n, mr_mip->modulus) == 0) return mr_mip->ndash;

    MR_IN(80)

    if (size(n) <= 2) {
        mr_berror(_MIPP_ MR_ERR_BAD_MODULUS);
        MR_OUT
        return (mr_small)0;
    }

    zero(mr_mip->w6);
    zero(mr_mip->w15);

    /* set a small negative QNR (on the assumption that n is prime!) */
    /* These defaults can be over-ridden                             */

    /* Did you know that for p=2 mod 3, -3 is a QNR? */

    mr_mip->pmod8 = remain(_MIPP_ n, 8);

    switch (mr_mip->pmod8) {
        case 0:
        case 1:
        case 2:
        case 4:
        case 6:
            mr_mip->qnr = 0; /* none defined */
            break;
        case 3:
            mr_mip->qnr = -1;
            break;
        case 5:
            mr_mip->qnr = -2;
            break;
        case 7:
            mr_mip->qnr = -1;
            break;
    }
    mr_mip->pmod9 = remain(_MIPP_ n, 9);

    mr_mip->NO_CARRY = FALSE;
    if (n->w[n->len - 1] >> M4 < 5) mr_mip->NO_CARRY = TRUE;

#ifdef MR_PENTIUM

    mr_mip->ACTIVE = FALSE;
    if (mr_mip->base != 0)
        if (MR_PENTIUM == n->len) mr_mip->ACTIVE = TRUE;
    if (MR_PENTIUM < 0) {
        if (n->len <= (-MR_PENTIUM)) mr_mip->ACTIVE = TRUE;
        if (logb2(_MIPP_ n) % mr_mip->lg2b == 0) mr_mip->ACTIVE = FALSE;
    }
#endif

#ifdef MR_DISABLE_MONTGOMERY
    mr_mip->MONTY = OFF;
#else
    mr_mip->MONTY = ON;
#endif

#ifdef MR_COMBA
    mr_mip->ACTIVE = FALSE;

    if (MR_COMBA == n->len && mr_mip->base == mr_mip->base2) {
        mr_mip->ACTIVE = TRUE;
#ifdef MR_SPECIAL
        mr_mip->MONTY = OFF; /* "special" modulus reduction */

#endif /* implemented in mrcomba.c    */
    }

#endif
    convert(_MIPP_ 1, mr_mip->one);
    if (!mr_mip->MONTY) { /* Montgomery arithmetic is turned off */
        copy(n, mr_mip->modulus);
        mr_mip->ndash = 0;
        MR_OUT
        return (mr_small)0;
    }

#ifdef MR_KCM

    /* test for base==0 & n->len=MR_KCM.2^x */

    mr_mip->ACTIVE = FALSE;
    if (mr_mip->base == 0) {
        nl = (int)n->len;
        while (nl >= MR_KCM) {
            if (nl == MR_KCM) {
                mr_mip->ACTIVE = TRUE;
                break;
            }
            if (nl % 2 != 0) break;
            nl /= 2;
        }
    }
    if (mr_mip->ACTIVE) {
        mr_mip->w6->len       = n->len + 1;
        mr_mip->w6->w[n->len] = 1;
        if (invmodp(_MIPP_ n, mr_mip->w6, mr_mip->w14) != 1) { /* problems */
            mr_berror(_MIPP_ MR_ERR_BAD_MODULUS);
            MR_OUT
            return (mr_small)0;
        }
    } else {
#endif
        mr_mip->w6->len   = 2;
        mr_mip->w6->w[0]  = 0;
        mr_mip->w6->w[1]  = 1; /* w6 = base */
        mr_mip->w15->len  = 1;
        mr_mip->w15->w[0] = n->w[0]; /* w15 = n mod base */
        if (invmodp(_MIPP_ mr_mip->w15, mr_mip->w6, mr_mip->w14) !=
            1) { /* problems */
            mr_berror(_MIPP_ MR_ERR_BAD_MODULUS);
            MR_OUT
            return (mr_small)0;
        }
#ifdef MR_KCM
    }
    copy(mr_mip->w14, mr_mip->big_ndash);
#endif

    mr_mip->ndash = mr_mip->base - mr_mip->w14->w[0]; /* = N' mod b */
    copy(n, mr_mip->modulus);
    mr_mip->check = OFF;
    mr_shift(_MIPP_ mr_mip->modulus, (int)mr_mip->modulus->len, mr_mip->pR);
    mr_mip->check = ON;
#ifdef MR_PENTIUM
    /* prime the FP stack */
    if (mr_mip->ACTIVE) {
        ndash = mr_mip->ndash;
        base  = mr_mip->base;
        magic *= base;
        ASM
        {
            finit fldcw WORD PTR control fld QWORD PTR ndash fld1 fld QWORD PTR
                base fdiv fld QWORD PTR magic
        }
    }
#endif
    nres(_MIPP_ mr_mip->one, mr_mip->one);
    MR_OUT

    return mr_mip->ndash;
}

void nres(_MIPD_ big x, big y)
{ /* convert x to n-residue format */
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif
    if (mr_mip->ERNUM) return;

    MR_IN(81)

    if (size(mr_mip->modulus) == 0) {
        mr_berror(_MIPP_ MR_ERR_NO_MODULUS);
        MR_OUT
        return;
    }
    copy(x, y);
    divide(_MIPP_ y, mr_mip->modulus, mr_mip->modulus);
    if (size(y) < 0) add(_MIPP_ y, mr_mip->modulus, y);
    if (!mr_mip->MONTY) {
        MR_OUT
        return;
    }
    mr_mip->check = OFF;

    mr_shift(_MIPP_ y, (int)mr_mip->modulus->len, mr_mip->w0);
    divide(_MIPP_ mr_mip->w0, mr_mip->modulus, mr_mip->modulus);
    mr_mip->check = ON;
    copy(mr_mip->w0, y);

    MR_OUT
}

void redc(_MIPD_ big x, big y)
{ /* Montgomery's REDC function p. 520 */
    /* also used to convert n-residues back to normal form */
    mr_small carry, delay_carry, m, ndash;  //,*w0g,*mg;

#ifdef MR_ITANIUM
    mr_small tm;
#endif
#ifdef MR_WIN64
    mr_small tm, tr;
#endif
    int i, j, rn, rn2;
    big w0, modulus;
#ifdef MR_NOASM
    union doubleword dble;
    mr_large         dbled;  //,ldres;
#endif
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif
    if (mr_mip->ERNUM) return;

    MR_IN(82)

    w0 = mr_mip->w0; /* get these into local variables (for inline assembly) */
    modulus = mr_mip->modulus;
    ndash   = mr_mip->ndash;

    copy(x, w0);
    if (!mr_mip->MONTY) {
        /*#ifdef MR_CELL
                mod256(_MIPP_ w0,w0);
        #else */
        divide(_MIPP_ w0, modulus, modulus);
        /* #endif */
        copy(w0, y);
        MR_OUT
        return;
    }
    delay_carry = 0;
    rn          = (int)modulus->len;
    rn2         = rn + rn;
#ifndef MR_SIMPLE_BASE
    if (mr_mip->base == 0) {
#endif
#ifndef MR_NOFULLWIDTH
        //   mg=modulus->w;
        //   w0g=w0->w;
        for (i = 0; i < rn; i++) {
            /* inline - substitutes for loop below */
#if INLINE_ASM == 1
            ASM cld ASM mov cx, rn ASM mov si, i ASM shl si,
                1
#ifdef MR_LMM
                ASM push ds ASM push es ASM les bx,
                DWORD PTR w0g ASM add bx, si ASM mov ax,
                es : [bx] ASM mul WORD PTR ndash ASM mov di,
                     ax ASM lds                          si,
                     DWORD PTR mg
#else
                ASM mov     bx,
                w0g ASM add bx,
                si ASM mov ax, [bx] ASM mul WORD PTR ndash ASM mov di,
                ax ASM mov si,
                mg
#endif
                             ASM push bp ASM xor
                         bp,
                     bp m1 : ASM lodsw ASM mul di ASM add ax,
                             bp ASM adc                   dx,
                             0
#ifdef MR_LMM
                             ASM add es : [bx],
                                          ax
#else
                        ASM add[bx],
                        ax
#endif
                                              ASM adc                     dx,
                                          0 ASM inc bx ASM inc bx ASM mov bp,
                                          dx ASM loop m1

                                              ASM pop bp ASM mov ax,
                                          delay_carry
#ifdef MR_LMM
                                              ASM add                 es : [bx],
                                                           ax ASM mov ax,
                                                           0 ASM adc  ax,
                                                           0 ASM add  es
                : [bx],
                  dx ASM pop es ASM pop ds
#else
                            ASM    add[bx],
                        ax ASM mov ax,
                        0 ASM adc  ax,
                        0 ASM      add[bx],
                        dx
#endif
                      ASM adc ax,
                  0 ASM mov   delay_carry,
                  ax
#endif
#if INLINE_ASM == 2
                      ASM cld ASM mov cx,
                  rn ASM mov          si,
                  i ASM shl           si,
                  2
#ifdef MR_LMM
                  ASM push ds ASM push es ASM les           bx,
                  DWORD PTR w0g ASM add                     bx,
                  si ASM mov                                eax,
                  es : [bx] ASM mul DWORD PTR ndash ASM mov edi,
                       eax ASM lds                          si,
                       DWORD PTR mg
#else
                        ASM mov                              bx,
                        w0g ASM add                          bx,
                        si ASM mov                           eax,
                        [bx] ASM mul DWORD PTR ndash ASM mov edi,
                        eax ASM mov                          si,
                        mg
#endif
                               ASM push ebp ASM xor
                           ebp,
                       ebp m1 : ASM lodsd ASM mul edi ASM add eax,
                                ebp ASM adc                   edx,
                                0
#ifdef MR_LMM
                                ASM add es : [bx],
                                             eax
#else
                                 ASM add[bx],
                                 eax
#endif
                                                 ASM adc edx,
                                             0 ASM add   bx,
                                             4 ASM mov   ebp,
                                             edx ASM loop m1

                                                 ASM pop ebp ASM mov eax,
                                             delay_carry
#ifdef MR_LMM
                                                 ASM add es : [bx],
                                                              eax ASM mov eax,
                                                              0 ASM adc   eax,
                                                              0 ASM add   es
                : [bx],
                  edx ASM pop es ASM pop ds
#else
                                     ASM     add[bx],
                                 eax ASM mov eax,
                                 0 ASM adc   eax,
                                 0 ASM       add[bx],
                                 edx
#endif
                      ASM adc eax,
                  0 ASM mov   delay_carry,
                  eax

#endif
#if INLINE_ASM == 3
                      ASM mov                           ecx,
                  rn ASM mov                            esi,
                  i ASM shl                             esi,
                  2 ASM mov                             ebx,
                  w0g ASM add                           ebx,
                  esi ASM mov                           eax,
                  [ebx] ASM mul DWORD PTR ndash ASM mov edi,
                  eax ASM mov                           esi,
                  mg ASM sub                            ebx,
                  esi ASM sub                           ebx,
                  4 ASM push ebp                        ASM xor ebp,
                  ebp                                   m1
                : ASM mov                               eax,
                  [esi] ASM add                         esi,
                  4 ASM mul edi ASM add                 eax,
                  ebp ASM mov                           ebp,
                  [esi + ebx] ASM adc                   edx,
                  0 ASM add                             ebp,
                  eax ASM adc                           edx,
                  0 ASM                                 mov[esi + ebx],
                  ebp ASM dec ecx ASM mov               ebp,
                  edx ASM jnz m1

                      ASM pop ebp ASM mov eax,
                  delay_carry ASM         add[esi + ebx + 4],
                  eax ASM mov             eax,
                  0 ASM adc               eax,
                  0 ASM                   add[esi + ebx + 4],
                  edx ASM adc             eax,
                  0 ASM mov               delay_carry,
                  eax

#endif
#if INLINE_ASM == 4
                      ASM("movl %0,%%ecx\n"
                          "movl %1,%%esi\n"
                          "shll $2,%%esi\n"
                          "movl %2,%%ebx\n"
                          "addl %%esi,%%ebx\n"
                          "movl (%%ebx),%%eax\n"
                          "mull %3\n"
                          "movl %%eax,%%edi\n"
                          "movl %4,%%esi\n"
                          "subl %%esi,%%ebx\n"
                          "subl $4,%%ebx\n"
                          "pushl %%ebp\n"
                          "xorl %%ebp,%%ebp\n"
                          "m1:\n"
                          "movl (%%esi),%%eax\n"
                          "addl $4,%%esi\n"
                          "mull %%edi\n"
                          "addl %%ebp,%%eax\n"
                          "movl (%%esi,%%ebx),%%ebp\n"
                          "adcl $0,%%edx\n"
                          "addl %%eax,%%ebp\n"
                          "adcl $0,%%edx\n"
                          "movl %%ebp,(%%esi,%%ebx)\n"
                          "decl %%ecx\n"
                          "movl %%edx,%%ebp\n"
                          "jnz m1\n"

                          "popl %%ebp\n"
                          "movl %5,%%eax\n"
                          "addl %%eax,4(%%esi,%%ebx)\n"
                          "movl $0,%%eax\n"
                          "adcl $0,%%eax\n"
                          "addl %%edx,4(%%esi,%%ebx)\n"
                          "adcl $0,%%eax\n"
                          "movl %%eax,%5\n"

                          :
                          : "m"(rn), "m"(i), "m"(w0g), "m"(ndash), "m"(mg),
                            "m"(delay_carry)
                          : "eax", "edi", "esi", "ebx", "ecx", "edx", "memory");
#endif

#ifndef INLINE_ASM
            /*        muldvd(w0->w[i],ndash,0,&m);    Note that after this time
             */
            m     = ndash * w0->w[i];
            carry = 0; /* around the loop, w0[i]=0    */

            for (j = 0; j < rn; j++) {
#ifdef MR_NOASM
                dble.d = (mr_large)m * modulus->w[j] + carry + w0->w[i + j];
                w0->w[i + j] = dble.h[MR_BOT];
                carry        = dble.h[MR_TOP];
#else
                muldvd2(m, modulus->w[j], &carry, &w0->w[i + j]);
#endif
            }
            w0->w[rn + i] += delay_carry;
            if (w0->w[rn + i] < delay_carry)
                delay_carry = 1;
            else
                delay_carry = 0;
            w0->w[rn + i] += carry;
            if (w0->w[rn + i] < carry) delay_carry = 1;
#endif
        }
#endif

#ifndef MR_SIMPLE_BASE
    } else
        for (i = 0; i < rn; i++) {
#ifdef MR_FP_ROUNDING
            imuldiv(w0->w[i], ndash, 0, mr_mip->base, mr_mip->inverse_base, &m);
#else
            muldiv(w0->w[i], ndash, 0, mr_mip->base, &m);
#endif
            carry = 0;
            for (j = 0; j < rn; j++) {
#ifdef MR_NOASM
                dbled = (mr_large)m * modulus->w[j] + carry + w0->w[i + j];
#ifdef MR_FP_ROUNDING
                carry = (mr_small)MR_LROUND(dbled * mr_mip->inverse_base);
#else
#ifndef MR_FP
                if (mr_mip->base == mr_mip->base2)
                    carry = (mr_small)(dbled >> mr_mip->lg2b);
                else
#endif
                    carry = (mr_small)MR_LROUND(dbled / mr_mip->base);
#endif
                w0->w[i + j] =
                    (mr_small)(dbled - (mr_large)carry * mr_mip->base);
#else
#ifdef MR_FP_ROUNDING
                carry =
                    imuldiv(modulus->w[j], m, w0->w[i + j] + carry,
                            mr_mip->base, mr_mip->inverse_base, &w0->w[i + j]);
#else
                carry = muldiv(modulus->w[j], m, w0->w[i + j] + carry,
                               mr_mip->base, &w0->w[i + j]);
#endif
#endif
            }
            w0->w[rn + i] += (delay_carry + carry);
            delay_carry = 0;
            if (w0->w[rn + i] >= mr_mip->base) {
                w0->w[rn + i] -= mr_mip->base;
                delay_carry = 1;
            }
        }
#endif
    w0->w[rn2] = delay_carry;
    w0->len    = rn2 + 1;
    mr_shift(_MIPP_ w0, (-rn), w0);
    mr_lzero(w0);

    if (mr_compare(w0, modulus) >= 0) mr_psub(_MIPP_ w0, modulus, w0);
    copy(w0, y);
    MR_OUT
}

/* "Complex" method for ZZn2 squaring */

void nres_complex(_MIPD_ big a, big b, big r, big i)
{
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif
    if (mr_mip->ERNUM) return;
    MR_IN(225)

    if (mr_mip->NO_CARRY &&
        mr_mip->qnr ==
            -1) { /* if modulus is small enough we can ignore carries, and use
                     simple addition and subtraction */
        /* recall that Montgomery reduction can cope as long as product is less
         * than pR */
#ifdef MR_COMBA
#ifdef MR_COUNT_OPS
        fpa += 3;
#endif
        if (mr_mip->ACTIVE) {
            comba_add(a, b, mr_mip->w1);
            comba_add(a, mr_mip->modulus, mr_mip->w2); /* a-b is p+a-b */
            comba_sub(mr_mip->w2, b, mr_mip->w2);
            comba_add(a, a, r);
        } else {
#endif
            mr_padd(_MIPP_ a, b, mr_mip->w1);
            mr_padd(_MIPP_ a, mr_mip->modulus, mr_mip->w2);
            mr_psub(_MIPP_ mr_mip->w2, b, mr_mip->w2);
            mr_padd(_MIPP_ a, a, r);
#ifdef MR_COMBA
        }
#endif
        nres_modmult(_MIPP_ r, b, i);
        nres_modmult(_MIPP_ mr_mip->w1, mr_mip->w2, r);
    } else {
        nres_modadd(_MIPP_ a, b, mr_mip->w1);
        nres_modsub(_MIPP_ a, b, mr_mip->w2);

        if (mr_mip->qnr == -2) nres_modsub(_MIPP_ mr_mip->w2, b, mr_mip->w2);

        nres_modmult(_MIPP_ a, b, i);
        nres_modmult(_MIPP_ mr_mip->w1, mr_mip->w2, r);

        if (mr_mip->qnr == -2) nres_modadd(_MIPP_ r, i, r);

        nres_modadd(_MIPP_ i, i, i);
    }
    MR_OUT
}

#ifndef MR_NO_LAZY_REDUCTION

/*

Lazy reduction technique for zzn2 multiplication - competitive if Reduction is
more expensive that Multiplication. This is true for pairing-based crypto. Note
that Lazy reduction can also be used with Karatsuba! Uses w1, w2, w5, and w6.

Reduction poly is X^2-D=0

(a0+a1.X).(b0+b1.X) = (a0.b0 + D.a1.b1) + (a1.b0+a0.b1).X

Karatsuba

   (a0.b0+D.a1.b1) + ((a0+a1)(b0+b1) - a0.b0 - a1.b1).X
*/

void nres_lazy(_MIPD_ big a0, big a1, big b0, big b1, big r, big i)
{
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif
    if (mr_mip->ERNUM) return;

    mr_mip->check = OFF;
#ifdef MR_COUNT_OPS
    fpc += 3;
    fpa += 5;
    if (mr_mip->qnr == -2) fpa++;
#endif

#ifdef MR_COMBA
    if (mr_mip->ACTIVE) {
        comba_mult(a0, b0, mr_mip->w0);
        comba_mult(a1, b1, mr_mip->w5);
    } else {
#endif
#ifdef MR_KCM
        if (mr_mip->ACTIVE) {
            kcm_mul(_MIPP_ a1, b1, mr_mip->w5); /* this destroys w0! */
            kcm_mul(_MIPP_ a0, b0, mr_mip->w0);
        } else {
#endif
            MR_IN(151)
            multiply(_MIPP_ a0, b0, mr_mip->w0);
            multiply(_MIPP_ a1, b1, mr_mip->w5);
#ifdef MR_COMBA
        }
#endif
#ifdef MR_KCM
    }
#endif

    if (mr_mip->NO_CARRY &&
        mr_mip->qnr ==
            -1) { /* if modulus is small enough we can ignore carries, and use
                     simple addition and subtraction */
#ifdef MR_COMBA
#ifdef MR_COUNT_OPS
        fpa += 2;
#endif
        if (mr_mip->ACTIVE) {
            comba_double_add(mr_mip->w0, mr_mip->w5, mr_mip->w6);
            comba_add(a0, a1, mr_mip->w1);
            comba_add(b0, b1, mr_mip->w2);
        } else {
#endif
            mr_padd(_MIPP_ mr_mip->w0, mr_mip->w5, mr_mip->w6);
            mr_padd(_MIPP_ a0, a1, mr_mip->w1);
            mr_padd(_MIPP_ b0, b1, mr_mip->w2);
#ifdef MR_COMBA
        }
#endif
    } else {
        nres_double_modadd(_MIPP_ mr_mip->w0, mr_mip->w5,
                           mr_mip->w6); /* w6 =  a0.b0+a1.b1 */
        if (mr_mip->qnr == -2)
            nres_double_modadd(_MIPP_ mr_mip->w5, mr_mip->w5, mr_mip->w5);
        nres_modadd(_MIPP_ a0, a1, mr_mip->w1);
        nres_modadd(_MIPP_ b0, b1, mr_mip->w2);
    }
    nres_double_modsub(_MIPP_ mr_mip->w0, mr_mip->w5,
                       mr_mip->w0); /* r = a0.b0+D.a1.b1 */

#ifdef MR_COMBA
    if (mr_mip->ACTIVE) {
        comba_redc(_MIPP_ mr_mip->w0, r);
        comba_mult(mr_mip->w1, mr_mip->w2, mr_mip->w0);
    } else {
#endif
#ifdef MR_KCM
        if (mr_mip->ACTIVE) {
            kcm_redc(_MIPP_ mr_mip->w0, r);
            kcm_mul(_MIPP_ mr_mip->w1, mr_mip->w2, mr_mip->w0);
        } else {
#endif
            redc(_MIPP_ mr_mip->w0, r);
            multiply(_MIPP_ mr_mip->w1, mr_mip->w2,
                     mr_mip->w0); /* w0=(a0+a1)*(b0+b1) */
#ifdef MR_COMBA
        }
#endif
#ifdef MR_KCM
    }
#endif

    if (mr_mip->NO_CARRY && mr_mip->qnr == -1) {
#ifdef MR_COMBA
        if (mr_mip->ACTIVE)
            comba_double_sub(mr_mip->w0, mr_mip->w6, mr_mip->w0);
        else
#endif
            mr_psub(_MIPP_ mr_mip->w0, mr_mip->w6, mr_mip->w0);
    } else
        nres_double_modsub(_MIPP_ mr_mip->w0, mr_mip->w6,
                           mr_mip->w0); /* (a0+a1)*(b0+b1) - w6 */

#ifdef MR_COMBA
    if (mr_mip->ACTIVE) {
        comba_redc(_MIPP_ mr_mip->w0, i);
    } else {
#endif
#ifdef MR_KCM
        if (mr_mip->ACTIVE) {
            kcm_redc(_MIPP_ mr_mip->w0, i);
        } else {
#endif
            redc(_MIPP_ mr_mip->w0, i);
            MR_OUT
#ifdef MR_COMBA
        }
#endif
#ifdef MR_KCM
    }
#endif

    mr_mip->check = ON;
}

#endif

#ifndef MR_STATIC

void nres_dotprod(_MIPD_ int n, big* x, big* y, big w)
{
    int i;
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif

    if (mr_mip->ERNUM) return;
    MR_IN(120)
    mr_mip->check = OFF;
    zero(mr_mip->w7);
    for (i = 0; i < n; i++) {
        multiply(_MIPP_ x[i], y[i], mr_mip->w0);
        mr_padd(_MIPP_ mr_mip->w7, mr_mip->w0, mr_mip->w7);
    }
    copy(mr_mip->pR, mr_mip->w6);
    /* w6 = p.R */
    divide(_MIPP_ mr_mip->w7, mr_mip->w6, mr_mip->w6);
    redc(_MIPP_ mr_mip->w7, w);

    mr_mip->check = ON;
    MR_OUT
}

#endif

void nres_negate(_MIPD_ big x, big w)
{
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif
    if (size(x) == 0) {
        zero(w);
        return;
    }
#ifdef MR_COMBA
    if (mr_mip->ACTIVE) {
        comba_negate(_MIPP_ x, w);
        return;
    } else {
#endif
        if (mr_mip->ERNUM) return;

        MR_IN(92)
        mr_psub(_MIPP_ mr_mip->modulus, x, w);
        MR_OUT

#ifdef MR_COMBA
    }
#endif
}

void nres_div2(_MIPD_ big x, big w)
{
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif

    MR_IN(198)
    copy(x, mr_mip->w1);
    if (remain(_MIPP_ mr_mip->w1, 2) != 0)
        add(_MIPP_ mr_mip->w1, mr_mip->modulus, mr_mip->w1);
    subdiv(_MIPP_ mr_mip->w1, 2, mr_mip->w1);
    copy(mr_mip->w1, w);

    MR_OUT
}

void nres_div3(_MIPD_ big x, big w)
{
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif

    MR_IN(199)
    copy(x, mr_mip->w1);
    while (remain(_MIPP_ mr_mip->w1, 3) != 0)
        add(_MIPP_ mr_mip->w1, mr_mip->modulus, mr_mip->w1);
    subdiv(_MIPP_ mr_mip->w1, 3, mr_mip->w1);
    copy(mr_mip->w1, w);

    MR_OUT
}

void nres_div5(_MIPD_ big x, big w)
{
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif

    MR_IN(208)
    copy(x, mr_mip->w1);
    while (remain(_MIPP_ mr_mip->w1, 5) != 0)
        add(_MIPP_ mr_mip->w1, mr_mip->modulus, mr_mip->w1);
    subdiv(_MIPP_ mr_mip->w1, 5, mr_mip->w1);
    copy(mr_mip->w1, w);

    MR_OUT
}

/* mod pR addition and subtraction */
#ifndef MR_NO_LAZY_REDUCTION

void nres_double_modadd(_MIPD_ big x, big y, big w)
{
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif
#ifdef MR_COMBA

    if (mr_mip->ACTIVE) {
        comba_double_modadd(_MIPP_ x, y, w);
        return;
    } else {
#endif

        if (mr_mip->ERNUM) return;
        MR_IN(153)

        mr_padd(_MIPP_ x, y, w);
        if (mr_compare(w, mr_mip->pR) >= 0) mr_psub(_MIPP_ w, mr_mip->pR, w);

        MR_OUT
#ifdef MR_COMBA
    }
#endif
}

void nres_double_modsub(_MIPD_ big x, big y, big w)
{
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif
#ifdef MR_COMBA

    if (mr_mip->ACTIVE) {
        comba_double_modsub(_MIPP_ x, y, w);
        return;
    } else {
#endif

        if (mr_mip->ERNUM) return;
        MR_IN(154)

        if (mr_compare(x, y) >= 0)
            mr_psub(_MIPP_ x, y, w);
        else {
            mr_psub(_MIPP_ y, x, w);
            mr_psub(_MIPP_ mr_mip->pR, w, w);
        }

        MR_OUT
#ifdef MR_COMBA
    }
#endif
}

#endif

void nres_modadd(_MIPD_ big x, big y, big w)
{ /* modular addition */
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif
#ifdef MR_COUNT_OPS
    fpa++;
#endif
#ifdef MR_COMBA

    if (mr_mip->ACTIVE) {
        comba_modadd(_MIPP_ x, y, w);
        return;
    } else {
#endif
        if (mr_mip->ERNUM) return;

        MR_IN(90)
        mr_padd(_MIPP_ x, y, w);
        if (mr_compare(w, mr_mip->modulus) >= 0)
            mr_psub(_MIPP_ w, mr_mip->modulus, w);

        MR_OUT
#ifdef MR_COMBA
    }
#endif
}

void nres_modsub(_MIPD_ big x, big y, big w)
{ /* modular subtraction */

#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif
#ifdef MR_COUNT_OPS
    fpa++;
#endif
#ifdef MR_COMBA
    if (mr_mip->ACTIVE) {
        comba_modsub(_MIPP_ x, y, w);
        return;
    } else {
#endif
        if (mr_mip->ERNUM) return;

        MR_IN(91)

        if (mr_compare(x, y) >= 0)
            mr_psub(_MIPP_ x, y, w);
        else {
            mr_psub(_MIPP_ y, x, w);
            mr_psub(_MIPP_ mr_mip->modulus, w, w);
        }

        MR_OUT
#ifdef MR_COMBA
    }
#endif
}

int nres_moddiv(_MIPD_ big x, big y, big w)
{ /* Modular division using n-residues w=x/y mod n */
    int gcd;
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif
    if (mr_mip->ERNUM) return 0;

    MR_IN(85)

    if (x == y) { /* Illegal parameter usage */
        mr_berror(_MIPP_ MR_ERR_BAD_PARAMETERS);
        MR_OUT

        return 0;
    }
    redc(_MIPP_ y, mr_mip->w6);
    gcd = invmodp(_MIPP_ mr_mip->w6, mr_mip->modulus, mr_mip->w6);

    if (gcd != 1)
        zero(w); /* fails silently and returns 0 */
    else {
        nres(_MIPP_ mr_mip->w6, mr_mip->w6);
        nres_modmult(_MIPP_ x, mr_mip->w6, w);
        /*    mad(_MIPP_ x,mr_mip->w6,x,mr_mip->modulus,mr_mip->modulus,w); */
    }
    MR_OUT
    return gcd;
}

void nres_premult(_MIPD_ big x, int k, big w)
{ /* multiply n-residue by small ordinary integer */
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif
    int sign = 0;
    if (k == 0) {
        zero(w);
        return;
    }
    if (k < 0) {
        k    = -k;
        sign = 1;
    }
    if (mr_mip->ERNUM) return;

    MR_IN(102)

    if (k <= 6) {
        switch (k) {
            case 1:
                copy(x, w);
                break;
            case 2:
                nres_modadd(_MIPP_ x, x, w);
                break;
            case 3:
                nres_modadd(_MIPP_ x, x, mr_mip->w0);
                nres_modadd(_MIPP_ x, mr_mip->w0, w);
                break;
            case 4:
                nres_modadd(_MIPP_ x, x, w);
                nres_modadd(_MIPP_ w, w, w);
                break;
            case 5:
                nres_modadd(_MIPP_ x, x, mr_mip->w0);
                nres_modadd(_MIPP_ mr_mip->w0, mr_mip->w0, mr_mip->w0);
                nres_modadd(_MIPP_ x, mr_mip->w0, w);
                break;
            case 6:
                nres_modadd(_MIPP_ x, x, w);
                nres_modadd(_MIPP_ w, w, mr_mip->w0);
                nres_modadd(_MIPP_ w, mr_mip->w0, w);
                break;
        }
        if (sign == 1) nres_negate(_MIPP_ w, w);
        MR_OUT
        return;
    }

    mr_pmul(_MIPP_ x, (mr_small)k, mr_mip->w0);
#ifdef MR_COMBA
#ifdef MR_SPECIAL
    comba_redc(_MIPP_ mr_mip->w0, w);
#else
    divide(_MIPP_ mr_mip->w0, mr_mip->modulus, mr_mip->modulus);
    copy(mr_mip->w0, w);
#endif
#else
    divide(_MIPP_ mr_mip->w0, mr_mip->modulus, mr_mip->modulus);
    copy(mr_mip->w0, w);
#endif

    if (sign == 1) nres_negate(_MIPP_ w, w);

    MR_OUT
}

void nres_modmult(_MIPD_ big x, big y, big w)
{ /* Modular multiplication using n-residues w=x*y mod n */
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif
    if ((x == NULL || x->len == 0) && x == w) return;
    if ((y == NULL || y->len == 0) && y == w) return;
    if (y == NULL || x == NULL || x->len == 0 || y->len == 0) {
        zero(w);
        return;
    }
#ifdef MR_COUNT_OPS
    fpc++;
#endif
#ifdef MR_COMBA
    if (mr_mip->ACTIVE) {
        if (x == y)
            comba_square(x, mr_mip->w0);
        else
            comba_mult(x, y, mr_mip->w0);
        comba_redc(_MIPP_ mr_mip->w0, w);
    } else {
#endif
#ifdef MR_KCM
        if (mr_mip->ACTIVE) {
            if (x == y)
                kcm_sqr(_MIPP_ x, mr_mip->w0);
            else
                kcm_mul(_MIPP_ x, y, mr_mip->w0);
            kcm_redc(_MIPP_ mr_mip->w0, w);
        } else {
#endif
#ifdef MR_PENTIUM
            if (mr_mip->ACTIVE) {
                if (x == y)
                    fastmodsquare(_MIPP_ x, w);
                else
                    fastmodmult(_MIPP_ x, y, w);
            } else {
#endif
                if (mr_mip->ERNUM) return;

                MR_IN(83)

                mr_mip->check = OFF;
                multiply(_MIPP_ x, y, mr_mip->w0);
                redc(_MIPP_ mr_mip->w0, w);
                mr_mip->check = ON;
                MR_OUT
#ifdef MR_COMBA
            }
#endif
#ifdef MR_KCM
        }
#endif
#ifdef MR_PENTIUM
    }
#endif
}

/* Montgomery's trick for finding multiple   *
 * simultaneous modular inverses             *
 * Based on the observation that             *
 *           1/x = yz*(1/xyz)                *
 *           1/y = xz*(1/xyz)                *
 *           1/z = xy*(1/xyz)                *
 * Why are all of Peter Montgomery's clever  *
 * algorithms always described as "tricks" ??*/

BOOL nres_double_inverse(_MIPD_ big x, big y, big w, big z)
{ /* find y=1/x mod n and z=1/w mod n */
    /* 1/x = w/xw, and 1/w = x/xw       */
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif
    MR_IN(145)

    nres_modmult(_MIPP_ x, w, mr_mip->w6); /* xw */

    if (size(mr_mip->w6) == 0) {
        mr_berror(_MIPP_ MR_ERR_DIV_BY_ZERO);
        MR_OUT
        return FALSE;
    }
    redc(_MIPP_ mr_mip->w6, mr_mip->w6);
    redc(_MIPP_ mr_mip->w6, mr_mip->w6);
    invmodp(_MIPP_ mr_mip->w6, mr_mip->modulus, mr_mip->w6);

    nres_modmult(_MIPP_ w, mr_mip->w6, mr_mip->w5);
    nres_modmult(_MIPP_ x, mr_mip->w6, z);
    copy(mr_mip->w5, y);

    MR_OUT
    return TRUE;
}

BOOL nres_multi_inverse(_MIPD_ int m, big* x, big* w)
{ /* find w[i]=1/x[i] mod n, for i=0 to m-1 *
   * x and w MUST be distinct               */
    int i;
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif
    if (m == 0) return TRUE;
    if (m < 0) return FALSE;
    MR_IN(118)

    if (x == w) {
        mr_berror(_MIPP_ MR_ERR_BAD_PARAMETERS);
        MR_OUT
        return FALSE;
    }

    if (m == 1) {
        copy(mr_mip->one, w[0]);
        nres_moddiv(_MIPP_ w[0], x[0], w[0]);
        MR_OUT
        return TRUE;
    }

    convert(_MIPP_ 1, w[0]);
    copy(x[0], w[1]);
    for (i = 2; i < m; i++) nres_modmult(_MIPP_ w[i - 1], x[i - 1], w[i]);

    nres_modmult(_MIPP_ w[m - 1], x[m - 1],
                 mr_mip->w6); /* y=x[0]*x[1]*x[2]....x[m-1] */
    if (size(mr_mip->w6) == 0) {
        mr_berror(_MIPP_ MR_ERR_DIV_BY_ZERO);
        MR_OUT
        return FALSE;
    }

    redc(_MIPP_ mr_mip->w6, mr_mip->w6);
    redc(_MIPP_ mr_mip->w6, mr_mip->w6);

    invmodp(_MIPP_ mr_mip->w6, mr_mip->modulus, mr_mip->w6);

    /* Now y=1/y */

    copy(x[m - 1], mr_mip->w5);
    nres_modmult(_MIPP_ w[m - 1], mr_mip->w6, w[m - 1]);

    for (i = m - 2;; i--) {
        if (i == 0) {
            nres_modmult(_MIPP_ mr_mip->w5, mr_mip->w6, w[0]);
            break;
        }
        nres_modmult(_MIPP_ w[i], mr_mip->w5, w[i]);
        nres_modmult(_MIPP_ w[i], mr_mip->w6, w[i]);
        nres_modmult(_MIPP_ mr_mip->w5, x[i], mr_mip->w5);
    }

    MR_OUT
    return TRUE;
}

/* initialise elliptic curve */

void ecurve_init(_MIPD_ big a, big b, big p, int type)
{ /* Initialize the active ecurve    *
   * Asize indicate size of A        *
   * Bsize indicate size of B        */
    int as;
#ifdef MR_OS_THREADS
    miracl* mr_mip = get_mip();
#endif
    if (mr_mip->ERNUM) return;

    MR_IN(93)

#ifndef MR_NO_SS
    mr_mip->SS = FALSE; /* no special support for super-singular curves */
#endif

    prepare_monty(_MIPP_ p);

    mr_mip->Asize = size(a);
    if (mr_abs(mr_mip->Asize) == MR_TOOBIG) {
        if (mr_mip->Asize >= 0) { /* big positive number - check it isn't minus
                                     something small */
            copy(a, mr_mip->w1);
            divide(_MIPP_ mr_mip->w1, p, p);
            subtract(_MIPP_ p, mr_mip->w1, mr_mip->w1);
            as = size(mr_mip->w1);
            if (as < MR_TOOBIG) mr_mip->Asize = -as;
        }
    }
    nres(_MIPP_ a, mr_mip->A);

    mr_mip->Bsize = size(b);
    if (mr_abs(mr_mip->Bsize) == MR_TOOBIG) {
        if (mr_mip->Bsize >= 0) { /* big positive number - check it isn't minus
                                     something small */
            copy(b, mr_mip->w1);
            divide(_MIPP_ mr_mip->w1, p, p);
            subtract(_MIPP_ p, mr_mip->w1, mr_mip->w1);
            as = size(mr_mip->w1);
            if (as < MR_TOOBIG) mr_mip->Bsize = -as;
        }
    }

    nres(_MIPP_ b, mr_mip->B);
#ifdef MR_EDWARDS
    mr_mip->coord = MR_PROJECTIVE; /* only type supported for Edwards curves */
#else
#ifndef MR_AFFINE_ONLY
    if (type == MR_BEST)
        mr_mip->coord = MR_PROJECTIVE;
    else
        mr_mip->coord = type;
#else
    if (type == MR_PROJECTIVE) mr_berror(_MIPP_ MR_ERR_NOT_SUPPORTED);
#endif
#endif
    MR_OUT
    return;
}
